### This is how to summarize the results from the MCMCs: the results are in Table 2 of the paper. 

dd<-c("C:/XunCao/Research/Portfolio/mcmc_net_polity/")
setwd(dd)

Year<-2001:2005 
i<-1

OUT<-read.table(paste(Year[i], "try 2/OUT", sep=""), header=T)


par(mfrow=c(4,4))                           #examine marginal mixing
for(i in c(4:19)) { plot(OUT[,i],type="l", main=colnames(OUT)[i]); abline(h=0) }

PS<-OUT[OUT$scan>round(max(OUT$scan)/2),-(1:3)]  #posterior samples, dropping 
                                                 #the first half of the chain 
                                                 #to allow for burn in

#gives mean, std dev, and .025,.5,.975 quantiles
M.SD.Q<-rbind( apply(PS,2,mean),apply(PS,2,sd),
                apply(PS,2,quantile,probs=c(.05,.5,.95)) )
print(M.SD.Q)


# distance<-latent<-corgrow<-corrup.s<-gdp.s<-gdppc.s<-corrup.r<-capopen.r<-gdp.r<-gdppc.r<-NULL
burn<-200  ## ignore the first 200 iterations. 


## for Xd: we have, in exact order, distance between capital cities (-), and correlation of growth rate (-) --- expected sign in the (). 
## sender effects: corruption (notice that higher values for less corrupted governments); GDP; GDP per caipta; 
## receiver effects: corruption; capital market opennes; GDP;  GDP per caipta;

for (i in 1:5){
out<-read.table(paste(Year[i], "try 2/OUT", sep=""), header=T)
Sum.post<-NULL
for (j in c(4:dim(out)[2])){   sum.post<-c(quantile(out[,j][-c(1:burn)], c(.025, .975))[1], 
                               mean(out[,j][-c(1:burn)]), 
                               quantile(out[,j][-c(1:burn)], c(.025, .975))[2], 
                               sd(out[,j][-c(1:burn)]))
                               Sum.post<-rbind(Sum.post, sum.post)
                               }
Sum.post<-round(Sum.post, digits=5)
colnames(Sum.post)<-c("2.5%", "Mean", "97.5%", "St.Dev")
rownames(Sum.post)<-colnames(out)[c(4:dim(out)[2])]
# Sum.post
dput(Sum.post, paste(Year[i], "posteriors", sep=""))
}

post.2001<-dget(paste("2001", "posteriors", sep=""))
post.2002<-dget(paste("2002", "posteriors", sep=""))
post.2003<-dget(paste("2003", "posteriors", sep=""))
post.2004<-dget(paste("2004", "posteriors", sep=""))
post.2005<-dget(paste("2005", "posteriors", sep=""))



### 90%
for (i in 1:5){
out<-read.table(paste(Year[i], "try 2/OUT", sep=""), header=T)
Sum.post<-NULL
for (j in c(4:dim(out)[2])){   sum.post<-c(quantile(out[,j][-c(1:burn)], c(.05, .95))[1], 
                               mean(out[,j][-c(1:burn)]), 
                               quantile(out[,j][-c(1:burn)], c(.05, .95))[2], 
                               sd(out[,j][-c(1:burn)]))
                               Sum.post<-rbind(Sum.post, sum.post)
                               }
Sum.post<-round(Sum.post, digits=5)
colnames(Sum.post)<-c("5%", "Mean", "95%", "St.Dev")
rownames(Sum.post)<-colnames(out)[c(4:dim(out)[2])]
# Sum.post
dput(Sum.post, paste(Year[i], "posteriors", sep=""))
}

post.2001<-dget(paste("2001", "posteriors", sep=""))
post.2002<-dget(paste("2002", "posteriors", sep=""))
post.2003<-dget(paste("2003", "posteriors", sep=""))
post.2004<-dget(paste("2004", "posteriors", sep=""))
post.2005<-dget(paste("2005", "posteriors", sep=""))


distance<-rbind(post.2001[1,1:3], post.2002[1,1:3], post.2003[1,1:3], post.2004[1,1:3],post.2005[1,1:3])
latent  <-rbind(rep(NA, 3),       post.2002[2,1:3], post.2003[2,1:3], post.2004[2,1:3],post.2005[2,1:3])
corgrow <-rbind(post.2001[2,1:3], post.2002[3,1:3], post.2003[3,1:3], post.2004[3,1:3],post.2005[3,1:3])
rho     <-rbind(post.2001[18,1:3], post.2002[19,1:3], post.2003[19,1:3], post.2004[19,1:3],post.2005[19,1:3])

pol.s   <-rbind(post.2001[4,1:3], post.2002[5,1:3], post.2003[5,1:3], post.2004[5,1:3],post.2005[5,1:3])
corrup.s<-rbind(post.2001[5,1:3], post.2002[6,1:3], post.2003[6,1:3], post.2004[6,1:3],post.2005[6,1:3])
capop.s <-rbind(post.2001[6,1:3], post.2002[7,1:3], post.2003[7,1:3], post.2004[7,1:3],post.2005[7,1:3])
gdp.s   <-rbind(post.2001[7,1:3], post.2002[8,1:3], post.2003[8,1:3], post.2004[8,1:3],post.2005[8,1:3])
gdppc.s <-rbind(post.2001[8,1:3], post.2002[9,1:3], post.2003[9,1:3], post.2004[9,1:3],post.2005[9,1:3])

pol.r   <-rbind(post.2001[9,1:3],  post.2002[10,1:3], post.2003[10,1:3], post.2004[10,1:3],post.2005[10,1:3])
corrup.r<-rbind(post.2001[10,1:3], post.2002[11,1:3], post.2003[11,1:3], post.2004[11,1:3],post.2005[11,1:3])
capop.r <-rbind(post.2001[11,1:3], post.2002[12,1:3], post.2003[12,1:3], post.2004[12,1:3],post.2005[12,1:3])
gdp.r   <-rbind(post.2001[12,1:3], post.2002[13,1:3], post.2003[13,1:3], post.2004[13,1:3],post.2005[13,1:3])
gdppc.r <-rbind(post.2001[13,1:3], post.2002[14,1:3], post.2003[14,1:3], post.2004[14,1:3],post.2005[14,1:3])


Post<-list(distance, latent, corgrow, rho,      
pol.s,
corrup.s,
capop.s,
gdp.s,   
gdppc.s, 
pol.r,
corrup.r,
capop.r, 
gdp.r,   
gdppc.r)

p.names<-c("Distance", "Lagged Latent Positions", "Correlation of GDP growth", expression("Reciprocity: "*rho),      
"Sender: polity",
"Sender: transparency",
"Sender: captial market openness",
"Sender: gdp",   
"Sender: gdp per capita", 
"Receiver: polity",
"Receiver: transparency",
"Receiver: captial market openness", 
"Receiver: gdp",   
"Receiver: gdp per capita")

par(mfrow=c(7,2))
par(mar=c(2,3,1.5,1.5))
par(cex=.7)
par(mgp=c(1.75,.75,0))

for (i in c(2, 4, 1, 3, 5, 10, 6, 11, 7, 12, 8, 13, 9, 14)){

daa<-Post[[i]]
daa.n<-p.names[i]

add<-.075*abs(min(daa[,1], na.rm=T)-max(daa[,3], na.rm=T))
plot(1:5, main="", xaxt="n", daa[,2], xlab="", type="n", ylab="", ylim=c(min(daa[,1], na.rm=T)-add, max(daa[,3], na.rm=T)+add)); abline(h=0, col="black", lwd=1)
title(main = list(daa.n, cex=1, font=2))
segments(1:5, daa[,1], 1:5, daa[,3], col="black", lwd=1.5); points(1:5, daa[,2], pch=19, cex=1)
axis(1, 1:5, as.character(Year[1:5]),font=1)

}


## sender effects of domestic institutions: more than polity this time. 
# par(mfrow=c(2,2))
par(mfrow=c(1,1))
par(mar=c(2,3,1.5,1.5))
par(cex=.7)
par(mgp=c(1.75,.75,0))

i<-10 ### so this is the receiver effect of domestic institutions: 

daa<-Post[[i]]
daa.n<-p.names[i]

add<-.075*abs(min(daa[,1], na.rm=T)-max(daa[,3], na.rm=T))
plot(1:5, main="", xaxt="n", daa[,2], xlab="", type="n", ylab="", ylim=c(min(daa[,1], na.rm=T)-add, max(daa[,3], na.rm=T)+add)); abline(h=0, col="black", lwd=1)
# title(main = list(daa.n, cex=1, font=2))
segments(1:5, daa[,1], 1:5, daa[,3], col="black", lwd=2); points(1:5, daa[,2], pch=19, cex=1.5)
axis(1, 1:5, as.character(Year[1:5]),font=1)




## dyadic only:

par(mfrow=c(2,2))
#par(mfrow=c(1,1))
par(mar=c(2,3,1.5,1.5))
par(cex=.7)
par(mgp=c(1.75,.75,0))

for (i in c(2, 4, 1, 3)){

daa<-Post[[i]]
daa.n<-p.names[i]

add<-.075*abs(min(daa[,1], na.rm=T)-max(daa[,3], na.rm=T))
plot(1:5, main="", xaxt="n", daa[,2], xlab="", type="n", ylab="", ylim=c(min(daa[,1], na.rm=T)-add, max(daa[,3], na.rm=T)+add)); abline(h=0, col="black", lwd=1)
title(main = list(daa.n, cex=1, font=2))
segments(1:5, daa[,1], 1:5, daa[,3], col="black", lwd=1.5); points(1:5, daa[,2], pch=19, cex=1)
axis(1, 1:5, as.character(Year[1:5]),font=1)

}









## sender and receiver: 
par(mfrow=c(5,2))
par(mar=c(2,3,1.5,1.5))
par(cex=.7)
par(mgp=c(1.75,.75,0))

for (i in c(5, 10, 6, 11, 7, 12, 8, 13, 9, 14)){

daa<-Post[[i]]
daa.n<-p.names[i]

add<-.075*abs(min(daa[,1], na.rm=T)-max(daa[,3], na.rm=T))
plot(1:5, main="", xaxt="n", daa[,2], xlab="", type="n", ylab="", ylim=c(min(daa[,1], na.rm=T)-add, max(daa[,3], na.rm=T)+add)); abline(h=0, col="black", lwd=1)
title(main = list(daa.n, cex=1, font=2))
segments(1:5, daa[,1], 1:5, daa[,3], col="black", lwd=1.6); points(1:5, daa[,2], pch=19, cex=1)
axis(1, 1:5, as.character(Year[1:5]),font=1)

}



library(xtable)

for (i in c(2)){
out<-read.table(paste(Year[i], "try 2/OUT", sep=""), header=T)
Sum.post<-NULL
for (j in c(4:dim(out)[2])){   sum.post<-c(quantile(out[,j][-c(1:burn)], c(.05, .95))[1], 
                               mean(out[,j][-c(1:burn)]), 
                               quantile(out[,j][-c(1:burn)], c(.05, .95))[2], 
                               sd(out[,j][-c(1:burn)]))
                               Sum.post<-rbind(Sum.post, sum.post)
                               }
Sum.post<-round(Sum.post, digits=5)
colnames(Sum.post)<-c("2.5%", "Mean", "97.5%", "St.Dev")
rownames(Sum.post)<-colnames(out)[c(4:dim(out)[2])]
# Sum.post
# dput(Sum.post, paste(Year[i], "posteriors", sep=""))
}
xtable(Sum.post[, -4], digits=c(3, 3, 3, 3))



######
par(mfrow=c(5,2))
par(mar=c(2,3,1.5,1.5))
par(cex=.7)
par(mgp=c(1.75,.75,0))



for (i in c(5, 10, 6, 11, 7, 12, 8, 13, 9, 14)){

daa<-Post[[i]]
daa.n<-p.names[i]

add<-.075*abs(min(daa[,1], na.rm=T)-max(daa[,3], na.rm=T))
plot(1:5, main="", xaxt="n", daa[,2], xlab="", type="n", ylab="", ylim=c(min(daa[,1], na.rm=T)-add, max(daa[,3], na.rm=T)+add)); abline(h=0, col="black", lwd=2)
title(main = list(daa.n, cex=1.1, font=2))
segments(1:5, daa[,1], 1:5, daa[,3], col="black", lwd=2); points(1:5, daa[,2], pch=19, cex=1.1)
axis(1, 1:5, as.character(Year[1:5]),font=2)

}










### 
